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Abstract. Using elementary methods, we define and derive a particular weighted average of the trapezoidal 
and composite trapezoidal rules and show that this approximation, as well as its composite, is straightforward in 

computation. This approximation and its composite, in their general forms, are shown to have predictable error 
patterns; thus, an extrapolation method can be used to increase the accuracy. We then derive the necessary 
weights to use an extrapolation method to reduce error and converge more quickly than Romberg integration by 
allowing for improved accuracy with fewer necessary subintervals. The procedure necessary to implement this 
alternative method is then carefully described, followed by two examples. 
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1. Introduction 

A variety of numerical integration methods exist with diverse rates of convergence, codability, and ease of compu- 
tation, and with varying requirements, both of the function to be integrated as well as the information required 

for approximation. In this paper, we will examine a new method of integration similar in application and scope to 
the Romberg method of integration, but with relatively faster rates of convergence because of fewer information 
requirements. 

After providing motivation for its construction, we introduce a form of the initial approximation, denoted A, 
with emphasis on the ease of this approximation's calculation. We derive A, a particular combination of the 

trapezoidal and composite trapezoidal rules, and then the remaining error of A, denoted a. This error is derived 
by weighting the errors of the trapezoidal rule and composite trapezoidal rule by the same weights used in the 
original approximation. 

Having A, we derive its composite, A^, and determine its error, after which it becomes apparent that an extrap- 
olation technique could be applied. We then derive the form of ft, the extrapolation factor necessary to eliminate 
the first error term between A and Am- We then find the general form of these possible extrapolation factors 
that are necessary to eliminate that first error term for all feasible combinations of A and A^,- We perform the 
first extrapolation and then we find the next few extrapolation factors necessary to eliminate further error terms. 

Afterwards, we provide in detail the procedure necessary to appropriately use the extrapolation factors and 
approximations to arrive at the most accurate approximation, before finally including two examples on this 
technique's implementation. 



2. Motivation 

Consider the following graph, Figure 1, of a monotonically increasing function. Often such a graph is used to 
introduce numerical integration concepts, particularly when the underestimate and overestimate rectangles are 
included, and the stepsizes, denoted h, are consistent as they arc here. As is often the case, we define the stepsize 
as h= {b — a)/n where & — a is the entire width of the integral, and n is the number of subintervals being used. 
Also, we define the sum of these underestimate rectangles as Ui, and the sum of these overestimate rectangles 
as Oi. 



2 



M. BRANDON YOUNGBERG 




Now consider if we were to segregate the differences between the underestimate and overestimate rectangles and 
to stack these differences. Because the stepsizes are equal, we know that the width of this stack will have the 
value h. Additionally, if we include the height of the lowest value for f{x), which here is at /(a), we have a 
column that is exactly the same height and width as the largest overestimate rectangle. This concept is depicted 
in Figure 2 below. 




1 1 1 — 5, 

Figure 2. 



Now consider if we were to divide the width of this new rectangle by the same number of underestimate or over- 
estimate rectangles that were in the original graph. We could then adjust the original function appropriately to 
fit in this new rectangle. The original values for f{x) would be in their original position, but now h = {b~a)/n'^. 
This concept is depicted in the following graph. Figure 3. 
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Figure 3. 



Intuitively, here the underestimate rectangles could be regarded as underestimates of the error of the original 
underestimates. Adding these underestimates to the original underestimates provides a more accurate estimate. 
Similarly, these new overestimates are overestimates of the error of the original overestimates. We must subtract 
these from the original overestimates to provide a more accurate estimate. Now consider if we were to repeat 
this process an infinite number of times, each iteration more narrow than the previous. We could then sum these 
iterations for a more accurate underestimate and overestimate. This sum for the underestimate will be referred 
to as Uoo, and for the overestimate, Ooo- Our approximation will be the average of Uoo and Ooo- Additionally, 
it will be shown in the next section that, strangely, for this approximation to be valid, it is not necessary for the 
function to be monotonic at all! 



3. Derivation of the Initial Approximation and Initial Error 
For i e {0, 1, . . . , n}, set Xi — a + ih. Define Ui to be the first underestimate, so that Ui — h f{^i)i a-nd Oi 

n n — 1 

to be the first overestimate, so that Oi — h J2 fi^i)- For future reference note that if S = ^ f{xi), then 

4=1 1=1 



C/i + Oi =/j(/(a) + 2E + /(fe)), 
C/i-Oi-/i(/(a)-/(6)). 

Let Uoo and O^o denote the limit of the underestimates and overestimates, respectively. Finally, set H = hf{a). 
We then have 



(1) C/oo = C/i + - -^H, 

n — 1 n — 1 

(2) Ooo = 0^-^ + -^H. 

n + 1 n + 1 



The average of (1) and (2) is then 
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(3) 



2 

Ui+Oi 



Ui 



Oi 



n — 1 n + 1 



n 



n + 1 n — 1 



H 



(n2 - + Oi) + (n + l)Ui - (n - l)Oi 



2n 



— 1 



rn2(C/i + Oi) + n(C/i + Oi) - 2nH 



/in2(/(a) + 2S + f{h)) + /in(/(a) - - 2hnf{a) 



;(/(a) + 2E + /(6) 



/in 



(/(«) + /(&)) 



n2- 1 



Now note that f (/(a) + 2S + /(6)) = ^{f{a) + 2S + is the estimate of /J /(a;) da; coming from the 

composite trapezoidal rule and that ^(/(a) + f{b)) = ^^{f{a) + f{b)) is the estimate of j'' f{x) dx coming 
from the ordinary trapezoidal rule. 

To find our error, we can use Taylor expansion to derive the trapezoidal rule. We expand around a and b, and 
then we combine the expressions to obtain 



Jjix) dx = J [/(a) + m] + ^[/'(a) - f'm + ^[/"(a) + f"{b)] + 0{h% 

where h = b — a. We can then follow a well-known procedure to further use Taylor expansion and remove odd 
powers of h, which gives us the Euler-Maclaurin summation formula shown here: 



(4) fix) dx = ^[/(a) + m] + ^[/'(a) - />)] - [/"'(«) - /">)] + 0{h'). 

Subtracting the trapezoidal rule from (4) results in the trapezoidal rule error. Denote the first term of this error 

as Eti; such that Exi = (o) ^ / (b)]- Further discussion of the derivation the trapezoidal rule, its error, 

and the Euler-Maclaurin summation formula can be found in [2, 4, 5]. 

Now when we construct a composite, intermediate points cancel. To demonstrate this, consider the construction 
of a composite trapezoidal rule with two parts: 



r fix) dx = |[/(a:i) + f{x2)] + hf{x2) + f{x^)] 

+ - /'(^2)] + ^[/'(^2) - /'(X3)] + 0{h% 

Notice that in the last two terms, the f{x2) terms cancel. This occurs for each set of derivatives as they are 
always expressed as differences using this derivation method. Additionally, similar to as how the trapezoidal 
rule uses the entire width of the integral as its stepsize, the composite trapezoidal rule uses the entire width of 
each subinterval as its stepsize. If w is the width of each of these subintervals, then in the composite trapezoidal 
rule, h = nw/n, but nw = b — a. So if in the trapezoidal rule error we replace h = {b — a) with h = {b — a)/n, 
we are left with the composite trapezoidal rule error. Then, this composite trapezoidal rule error is written as 



Denote the first term of this error as Eci- Then, notice that Eqi, when multiplied by n^, is exactly the same as 
Eti- Then we can rewrite expression (3) to eliminate these first error terms and to derive the remaining error, 
denoted as a, as follows: 



ALTERNATIVE TO THE ROMBERG METHOD OF ESTIMATING THE DEFINITE INTEGRAL 



5 



h 



(/(a) + 2I] + /(6) 



hn 



f{x) dx + Ec 



n 


2-1 


Ec 






- 1 



f{x) dx + Et 



n 



f m dx + ^^^^[/"' (a) - f"'{b)] + Oih') 

J a 



{h-af 
6!n4 



2(3!)n2 



[/'(«)-/'(&)] 



■n? - 1 



+ - 



2(3!)n^ 



' f{x) dx + ^^^Ifia) - f"'{b)] + 0{{b - af) 



in' - 1) 



J a 



dx 



-n2(n2-l) 



6!n4 



V? - 1 

Then, A = f f{x) dx - ^^^^[f {a) - /" (b)] + O {g{h)), and 



(5) 

where 



0(5(/i)). 



= c 



2/^6 _ (5 _ a)e 



-n2(6-a)6-(6-a)« 



and where c is the arbitrary constant. It will be shown in the following section that g{h) will be able to be further 
simplified quite nicely, and it has been purposefully left in this form to more easily facilitate this simplification. 

4. Derivation of the Composite and the Extrapolation Factors 

Denote the composite approximation as A^, where m is the number of subintervals used in each part of the 
composite. We can clarify the relationship between the original approximation's error, A, and that of the 
composite approximation. Am, by first observing that the original approximation's error is such that 

a = - ^ttf [/"' («) - /"' (b)] + O {g{hj) , 



6!n2 

n^(b-a)^ 
6!n4 



6! 



[/ (a)-/ {b)]+0{g{h))., 

[/"'(«)-/"'(&)] +0 (.9 W). 



Consider, again, if we were to construct a composite with k parts of equal width w. Then, since each w is 
equivalent, kw = b — a and km = n. This allows us to now define our stepsize for the entire integral as 

h = w/m = kw/km ~ (b — a)/n. Then, each part of our composite would have the following error, denoted cKjt^, 
where the subscript m corresponds to the number of subintervals used in each part of the composite: 



6!m^ 



■[/ (a)-/ {b)]+0{j{h)), 



where 



V? — 1 



= c 



-m'^{b-af -{b-af 



Now the composite error for the entire integral becomes 
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(6) = ^g,^/^ [/"{a) - f"\b)] + O m) . 

It is now apparent that an extrapolation is possible since we will have two approximations with comparable 
errors, the original and a composite. Then we must first solve for the first weight, or extrapolation factor, used 
in our extrapolation process. We will use Q to denote these extrapolation factors, with subscripts to carefully 
track for which extrapolation each O is used. Calculation of fln,mj where the subscript n denotes the total number 
of subintervals in the entire integral and where m denotes the number of subintervals used in the construction 
of the composite, is as follows: 



{b-aY 



if (a)-/ m=nr. 



w?{h-aY 
6!n4 



[/ («)-/ {b)] 



(7) 



n y 



6!n2 

" 
1 

This holds providing m divides n evenly and 2 < m < n/2, as a larger m disallows construction of composites 

with parts of equal numbers of subintervals. From this point forward, we shall denote A as and its error as 
a as a„, because we can regard A^ as a composite constructed of one part with n subintervals. 

We can now red^ice the error from 0{n'^h^) to 0{n^h^) provided we use A„ and any feasible A„i; providing 
the composite is constructed with parts of equal numbers of subintervals. This straightforward procedure is 
conducted as follows: 



{-)' 



f{x) dx — a. 



m 



/ f{x) dx - a„ 

J a 



1 



{-)' 



fix) dx + ^ l;^/^ [/■■ (a) - /■■ (6)] + O im) 



' fix) dx + ^^-^ [/"' (a) - /"' (6)] + O igih)) 
bin'' 



n Y 
mJ 



f 

J a 

I 

J a 



fix) dx + 



(-Ym-gih) 



= / fix) dx + 0{c 



i-m^rr - ■n?)ih - a)^ + (m^n^ + m^)ih - a) 



n'^iv? — m2) 



f 



fix) dx + Oin^K'). 



Additionally, note that further error terms are described generally by 0('n?h'^) where p is even. This is impor- 
tant to recognize, as an increase in the number of points available will not necessarily increase the accuracy by 
an expected amount. For instance, should the first error term be simply 0(/i^), then doubling the number of 
points will increase the accuracy by a factor of G4. However, doubling the number of available points when the 
first error term is described as 0(n^/i^) will only initially increase the accuracy by a factor of 16. Incidentally, 
doubling the number of available points would also allow for the use of an additional composite, which would 
ultimately increase the accuracy, as will be demonstrated in the following sections. 
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At this point, it is interesting to note that, after some algebra, we can rewrite our approximation An in the 
following form: 



(8) = ^ 



nfja) ^ 2n^E ^ nfjb) 
n + 1 n^ — 1 n+1 



Then, when n = 2, (8) is exactly Simpson's rule with (5) and (6) exactly Simpson's Rule error and composite 
Simpson's Rule error, respectively. Additionally, when n = 3, the same is true concerning Simpson's 3/8 Rule 
and errors as well. However, beyond the case where n = 3, (8) ceases to follow further Newton-Cotes formulas 
exactly. Further discussion of the Newton-Cotes formulas and errors, including Simpson's Rule and Simpson's 
3/8 Rule can be found in [1, 6]. 

The extrapolation factors necessary to remove the remaining error terms are very simple to compute. Consider 

two separate composites, one comprised of parts of x subintervals. denoted A^, and one comprised of parts of y 
subintervals, denoted Ay. Then let the remaining error terms of these composites to be and ay, respectively. 
As shown in (7), we will have two separate extrapolation factors, and Then we know that 



and also that 



This allows us to then compute 

(9) ay = = (l) a^. 

We can now extend the use of the extrapolation factors' subscripts to denote combinations of previously-used 
extrapolation factors from which they result, such as: 



^- = ^ = (-)'(-)'=(-)'• 



Now consider if we had a third composite composed of parts of z subintervals. The associated error would then 
be denoted a^- We could then calculate a^, similarly as to in (9), giving us 

Q ^n,y f ^\ 

a^ — ^^z,y^y — ^y — 1 I 



Furthermore, 



Then generally. 



a. 



"•^ = ( y) ^V= (y) ) = ^z,y^y,xOiy = (^) a^ = 0^,^a^. 



)2 / 2 C ^ / Z ^ 



Therefore, any a can be combined with an appropriate extrapolation factor, fl, where the ft used is simply a 
combination of previously- used Q, in such a way as to equate this a with a^. If we allow z = n. where n is 
the total number of subintervals available in the entire integral, it becomes clear that we have, therefore, found 
the general form of all of the extrapolation factors necessary to reduce the error to the smallest value that the 
number of available composites will allow. Additionally, the general extrapolation factor Q.z,a is sufficient to 
convert all remaining error terms of our new approximation from Oia^h'') to 0{z'^h''). 
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5. Procedure 

Now that we have specified how to compute An, Am, and the extrapolation factor O, we use the following pro- 
cedure to produce the most accurate estimate available. Throughout this procedure, we must keep careful track 
of the subscripts of our extrapolation factors, as well as how we apply these factors in the extrapolation. Failure 
to carefully note the number of subintervals used in the construction of a composite, and how its associated 
extrapolation factor is used, may cause error in calculation of further extrapolations. Table 1 at the end of this 
section is provided as an aid for following the procedure. 

The first step is to compute A^, using the average of (1) and (2) as described in section 3, as well as A^ 
for all available composites. The composites must be constructed with parts comprised of an equal number of 

subintervals. Though these composites, as well as the initial approximation, will likely be different, they all have 
their first error term of 0{m?h'^), where m is the number of subintervals used in each composite. 

Suppose there are four separate available composites, which we will denote as Aa, Ah, Ac, and A^, where the 
number of subintervals used in each part of each composite is a, b, c, and d respectively. As previously shown in 
(7), the O for the extrapolation between An and Aa will be {n/a)'^, between Aa and A^ will be (a/b)'^, between 
Ah and Ac will be {b/cf', and between Ac and A^ will be {c/d}^. 

Perform the first extrapolation for each composite. Remember that care must be used to track the extrapolation 

factors for use in further extrapolations and computations of Q. The numerator of the extrapolation factor corre- 
sponds to the approximation being subtracted in the extrapolation process. For instance, the first extrapolation 
will be conducted using the original approximation, A„ and a composite as follows: 




Notice again that the numerator of the extrapolation factor, in this case n, corresponds with the approximation 
being subtracted. The next approximation would be 



Aa,b ~ ^ 



^a,b 



1 



Continuing in this manner, we should have several new approximations, one for each extrapolation we have per- 
formed so far. Denote these new approximations An^a, ^o,6) ^6,0 and Ac^d- Now we compute new extrapolation 
factors. As was previously shown in (9), the extrapolation factor to reduce the error using approximations An^a 
and Aa,b is fin, 6 = {^/o-YiO'/b)'^ = (j^/^)^- Now perform the next set of extrapolations as shown here with An^a 
and Aafi'- 




Again, we have new approximations, An.h, Aa.c, and Ah.d- Compute the next extrapolation factors, ftn,c and 
i^a,d- Then compute the next set of extrapolations and obtain the next two approximations. This procedure 
conducted with approximations An^b and Aa,c using 17„^c = {n/a)^{a/b)'^{b/c)^ = {n/c)"^ is done so in the 
following manner: 



O A — A I. (~] Aa,c — Anfi 



^n,c -i 



2 



1 



Two new approximations result, An^c and Aa,d- Calculate the next extrapolation factor, Q.n,d = {n/a)'^{a/b)'^{b/c)^{c/d)'^ = 
{n/df. Then, 
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O , 4 , _ 4 (-,} Aa,d — ^n,c 

This provides the highest accuracy estimation we are able to achieve given the original approximation, An and 
four composites. Note that this final approximation will have its first error term as 0{'n?h}'^), because all error 
terms of the form 0{m?h^'^) have been replaced with the error corresponding to the composite composed of 
parts of n subintervals, the original approximation. 

The following table helps to keep track of the extrapolation procedure and to clarify the manner in which this 
procedure causes convergence. Additionally, when implementing this procedure it is tempting to further reduce 
the values for Q, or to express them in decimal form, but we have left them reduced only to the level used in this 
procedure, as doing so may help to reduce confusion when constructing further O. 



Table 1 

Error 0(n''h'^) 0(n^h'^) 0(n?h^) ' 0(n2/i«) Oiri'h^^) 





- 1 




-A 




- 1 




-A,, 




- 1 




-Ac 





An,b 




- 1 








iff- 


- 1 



5)- ■ ■ G) 



- 1 



- 1 



6. Examples 

Example 1: 

This example was chosen because we know with the use of extrapolation we can arrive at an error term that 
is a multiple of derivatives that have the same constant value. This will ensure that our solution is exact. 
Additionally, h = 1 which simplifies some of our calculations. 

Consider f{x) = x'' ~ 2x + 10. Then f{x) dx = 12500000. Assume, however, that we do not know the 
function, but are instead given eleven equally spaced points from x = to 10 and their corresponding values for 



X f{x) 





10 


1 


9 


2 


134 


3 


2191 


4 


16386 


5 


78125 


6 


279934 


7 


823539 


8 


2097146 


9 


4782961 


10 


9999990 
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Our first step is to calculate Aiq by first calculating Ucc and Ooc- However, note that b — a = 10 = n, 
and therefore, h = I. Also, note that H = hf{a) = 10. Then, as stated previously, we calculate Ux as the 
underestimate, which is h multiplied by the sum of the first 10 values of f{x), which gives us Ui = 8080435. 
Likewise, the overestimate becomes Oi = 18080415. Then we have 



(10) U^ = U, + ^ "^H = 8080435 + - • 10 = 8978250, 

n— In— 1 9 9 

(11) 0^ = 0,-^ + -^H = 18080415 - + ^ • 10 = 16436750. 

n + 1 n + 1 11 11 

Averaging (10) and (11) then yields the initial approximation, Aio = 12707500. 

Now we can construct Am- We have two options when n = 10, a composite constructed using two parts of five 
subintervals and a composite constructed using five parts of two subintervals. 

Though the composite using five parts would provide a more accurate estimate, for the purposes of this example 
we will start with the composite using two parts. In this case note that, for each of these parts, m = 5. This 
does not change the value of h, however, as the width of each part of the composite is also 5. It does, however, 
change the value of H. Because this composite is found by simply adding the results of the two separate approx- 
imations using five subintervals each, we will have a different /(a) for each of these parts, and thus a separate 
H. In the first part. Hi = hf{0) = 10. However, in the second part, H2 = hf{5) = 78125. We will also have 
two separate sets of overestimates and underestimates, one set for each part of which the composite is comprised. 

The first underestimate is found using values /(O) through /(4), and our first overestimate is found using values 

/(I) through /(5). Then our second underestimate is found using values /(5) through /(9), and our second 
overestimate is found by using values /(6) through /(lO). Then we have 



(12) U^ = Ui + - -^Hi = 18730 + - 7 • 10 = 23400, 

n— In— 1 4 4 

(13) Ooo = Oi - ^ + -^Hi = 96845 - ^ + ^ 10 = 80712.5. 

n+ln+1 66 

Averaging (12) and (13) gives us 52056.25, which is our approximation of the area of the first half of the entire 
integral. 



(14) [/oo = f/2 + - = 8061705 + ^^^^ _ ^ . 78125 = 9979475, 

n— In— 1 44 

(15) = O2 - + -^H2 = 17983570 - + 5 . 78125 = 15051412.5. 

n+ln+l 6 6 

Averaging (14) and (15) gives us 12515443.75, which is our approximation of the area of the second half of the 
entire integral. Adding this to the approximation of the first half gives us A^ = 12567500. 

Now that we have our initial approximation and a composite approximation, we can find the appropriate ex- 
trapolation factor r^io.5- As this is our first extrapolation, we know that the numerator of Oio,5 is the total 
number of subintervals in the original approximation and that the denominator is the number of subintervals 
used in each part of the composite. Therefore, in this case, f2io,5 = (n/m)^ — (10/5)^. 

We calculate our new approximation, ^10,5, as follows: 



— j Ag-^io ( — ] 12567500- 12707500 

lio 5 = — o = — o = 12520833.3. 

'10\ . fW 



57 ' V5 

Calculating a composite using parts containing two subintervals each gives us A2 = 12511500. We then find 
^5,2 = 12500833.3 using ^5,2 = riio,2/f^io,5 = (10/2)V(10/5)2 = (5/2)2. -yy^ ^^^e to further combine 

these two approximations, ^10,5 and ^5^2, in such a manner as to further reduce the error. As shown in (9), in 
this second extrapolation, fiio,2 is simply a combination of previous extrapolation factors such that: 
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^10,2 = f^lO, 5^^5,2 = 

Therefore, our next extrapolation will be ^10,2 as computed below: 

12500833.3- 12520833.3 

Aio 2 = — n = 12500000. 

'10' 

2" 

This estimate is exact because we have reduced the error to 0{nPh^). This error term contains the difference of 

the seventh derivative of the function evaluated at a and the seventh derivative evaluated at b. However, that 
derivative is just a constant, allowing f^'^\a) = /'•^•'(&), which then allows this error term to negate itself. 



Example 2: 

This example was chosen because it is continuously differentiable. and thus, using the alternative integration 
method presented in this paper, we will always have some non-zero error term. It also allows for four separate 
composites to be computed, and we can more readily observe the role of h, because in this case, h^l. 

Consider f{x) = sin{x). Then J^^ f{x) dx = —2. Assume, however, that we do not know the function, but are 
instead given 13 equally spaced points from a; = tt to 27r and their corresponding values for f{x), rounded to 10 
decimal points. 



X 


fix) 


TT 





137r/12 


-0.2588190451 


77r/6 


-0.5 


157r/12 


-0.7071067812 


47r/3 


-0.8660254038 


177r/12 


-0.9659258263 


37r/2 


1 


1977/12 


-0.9659258263 


57r/3 


-0.8660254038 


2l7r/12 


-0.7071067812 


ll7r/6 


-0.5 


237r/12 


-0.2588190451 


27r 






Again, we calculate A12 by first calculating U^x, and O^o- However, note that n = 12. and therefore, h = 
{b — a)/n = 7r/12. Also, note that H = hf{a) = 0. Then, as stated previously, we calculate Ux as the 
underestimate, which is h multiplied by the sum of the first 12 values of f{x), which gives us Ui = —1.9885637766. 
Likewise, the overestimate becomes Oi = —1.9885637766. Then we have 

, , -1.9885637766 12 

(16) Uoo = -1.9885637766 -t- — n ' ° " -2.1693423017, 

, , ^ -1.9885637766 12 

(17) 0<x> = -1.9885637766 -h — • = -1.8355973322. 

\ J 00 13 13 

Averaging (16) and (17) then yields the initial approximation, A12 = —2.0024698170. 

Now we can construct the various available composites. We have four options when n = 12, a composite con- 
structed using two, three, four, or six parts of six, four, three, or two subintervals each, respectively. Denoting 
each composite with the number of subintervals used in each part during its construction gives us A2, A3, A4, 
and Aq. We would expect A2 to be the most accurate but most difiicult to compute, as it is a composite of 
many parts, and Aq to be the least accurate but easiest to compute. 
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Following the same procedure as in example 1, is calculated as follows: 

The first underestimate is found using values f(7r) through f(177r/12), and our first overestimate is found us- 
ing values f(137r/12) through f(37r/2). Then our second underestimate is found using values f(37r/2) through 
f(237r/12), and our second overestimate is found by using values f(197r/12) through f(27r). Then we have 

, , -0.8633821944 6 

(18) Uoo = -0.8633821944 + = -1.0360586333, 

5 5 

, , -1.1251815822 6 

(19) Ooo = -1.1251815822 + ;^ • = -0.9644413562. 

Averaging (18) and (19) gives us —1.0002499947, which is our approximation of the area of the first half of the 
entire integral. Note that, as in the previous example, we must now calculate a new H, as f{a) in this second 
half differs from that in the first half. Then we have i?2 = /i/(37r/2) = -0.2617993878, and 

, , -1.1251815822 6 

(20) Uoo = -1.1251815822 + • -0.2617993878 = -1.0360586333, 

O 

, , ^ -0.8633821944 6 

(21) = -0.8633821944 + - • -0.2617993878 = -0.9644413562. 

Averaging (20) and (21) gives us —1.0002499947, which is our approximation of the area of the second half of 
the entire integral. Adding this to the first half approximation gives us Aq = —2.0004999894. 

Now that we have our initial approximation and a composite approximation, we can find the appropriate extrap- 
olation factor, 012,6- Again, as this is our first extrapolation, wc know that the numerator of ^12,6 is the total 
number of subintervals in the original approximation and that the denominator is the number of subintervals 
used in each iteration of the composite. Therefore, in this case, Oi2,6 = (n/m)^ = (12/6)^. 



We calculate our new approximation, Ai2fi, as follows: 



12V . . /12^ ^ 



^6-^12 ( Y 1 (-2.0004999894) - (-2.0024698170) 

I12 6 = ^ — o = — 9 = -1.9998433802. 

/12\ , /12' 



Using the same procedure to calculate the other available composites gives us A2 = —2.0000526243, ^3 = 
—2.0001193864, and A4 = —2.0002147374. Now that we have our composite estimates, all that is required is to 
follow the previously described method to calculate our fl and our extrapolations, and we can construct Table 
2 similar to how Table 1 was constructed in the procedure section. 



Error 

A12 

Ae 

A2 

A3 

Ai 



0{n''h^) 
-2.0024698170 
-2.0004999894 
-2.0000526243 
-2.0001193864 
-2.0002147374 



-1.9998433802 
-1.9999967037 
-1.9999992147 
-1.9999967923 



Table 2. 

-2.0000010844 
-2.0000000517 
-2.0000000221 



-1.9999999828 
-1.9999999985 



-2.0000000005 



Note that in both examples, /(a) was not the lowest value for f{x) over the entire integral. The terms "underes- 
timate" and "overestimate" are, in some cases, misnomers. If the function is strictly decreasing over the integral, 
nothing has to be changed in the calculation of An or Am- This would result in the underestimate being larger 
than the overestimate, but this should not be alarming because the results of this procedure are not affected. 

7. Comparison with Romberg Integration 

Let us consider the results of Example 1. If we had decided instead to use the Romberg method of integration 
and were unable to increase the number of points, we would have had an initial estimate of 50,000,000 using 
only the endpoints and the trapezoidal rule. Our second estimate would have used the center point, essentially 
constructing a composite trapezoidal rule of two parts, but would not have fared much better at 25,390,625. 
If we had then used the extrapolation method inherent in Romberg integration, our estimate would have been 
reduced to 17,187,500, but we would be unable to continue further at this point, with our best estimate having 
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an error of almost 4.7 million. 

This is perhaps not a fair assessment of the Romberg method, because the polynomial used is of high degree, 
which can result in large values for the 0{h'^) error term. Furthermore, the Romberg method is crippled by being 
unable to use the optimal number of points. Then let's allow the number of points to increase to 17 such that 
h = 10/16. Counting the original trapezoidal estimation, Romberg integration will be able to use five separate 
approximations. In this case, Romberg integration will be able to easily reach 0{h^) error to arrive at an exact 
answer. 

In fact, the Romberg method would also be able to arrive at this accuracy with only 9 points, or 8 subintcrvals. 
However, with those same 8 subintervals, the alternative integration method presented in this paper would also 
arrive at the exact value. Unfortunately though, Romberg integration cannot reach this same accuracy with any 
number of subintcrvals less than 8, or between 8 and 16. It is only able to effectively use composites constructed 
of parts where the number of subintervals used in each part are powers of 2. Considering only the cases where 
we are given 16 or fewer subintervals, our alternative method would be able to reach an accuracy of 0{'n?h^). 
In this example, we would arrive at the exact answer, given 6, 8, 10, 12, 14, 15 or 16 subintervals, with 12 and 
16 subintervals providing even higher degrees of accuracy. 

Let us then consider Example 2. With the given 13 points, or 12 subintervals. Romberg integration will only 
be able to arrive at an estimate of -1.998570731824, which is a result of 0{h^) error. It would make use of the 
approximations from the normal trapezoidal rule, the composite trapezoidal rule composed of two parts, and the 
composite trapezoidal rule composed of four parts. Any further composites that could be used are not available. 

Suppose, again, that this is not a fair assessment. Let's provide the Romberg method with an ideal number of 

subintervals. 32. In this case, the Romberg method will be able to reach the error term 0{h^^). As is commonly 
done when implementing Romberg integration, we could build the following table of results. 



7r/2 

7r/4 

7r/8 

7r/16 

7r/32 





-1.57079632679490 
-1.89611889793705 
-1.97423160194556 
-1.99357034377234 
-1.99839336097015 



om 

-2.09439510239320 
-2.00455975498443 
-2.00026916994839 
-2.00001659104794 
-2.00000103336942 



Table 3. 



-1.99857073182384 
-1.99998313094599 
-1.99999975245458 
-1.99999999619085 



0{h^) 

-2.00000554997968 
-2.00000001628805 
-2.00000000005968 



-1.99999999458730 
-1.99999999999604 



0(/(") 
-2.00000000000133 



The alternative method introduced in this paper can use the same composites, except for the one using 32 subin- 
tervals, as that composite is incorporated in the original approximation, A^2- However, because the Romberg 
method's initial estimations have error of 0{h^), and the method introduced in this paper has initial estimations 
with an error of O^ri^h^), these two methods are easily comparable. Similarly as to how Table 3 was created, 
using the integration method presented in this paper, we can create Table 4, which, in practice, may be useful 
in keeping track of the various values of our extrapolation factors and approximations. 



Table 4. 



Error 


0(n2/i4) 


0(r),2/(,'') 


0(r|2/t») 




0(n2fi,i2) 


A32 


-2.00034682466611 


-1.99999967732512 


-2.00000000125221 


-1.99999999998018 


-2.00000000000133 


A2 


-2.00000103336942 


-1.99999999619085 


-2.00000000005968 


-1.99999999999604 




Ai 


-2.00000414490512 


-1.99999993815840 


-2.00000000406904 






As 


-2.00001676514528 


-1.99999894949885 








A16 


-2.00007021208456 











Though the final approximations reached by both methods is exactly the same, this does not appear to always 
be the case in every similar situation. However, when the estimation methods provide different results using 
the same composites, it is currently unclear as to what degree the difi^erences, often at the 15th or further digit, 
could be due to round-off error. Additionally, we were able to reach this final approximation as the 15th ap- 
proximation, where the Romberg method uses 21 approximations to arrive at the same value. 

Furthermore, note that if we had any number of subintervals less than 32, we would not be able to achieve 
this level of accuracy with Romberg integration. In fact, we would not be able to achieve this level of accuracy 
again until we were given 64 subintervals. However, our alternative method would be able to reach this level of 
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accuracy given 12, 18, 20, 24, 28, or 30 subintervals, for those less than 32, and even an accuracy of 0{n'^h^^) 
given either 24 or 30 subintervals. A thorough discussion of the principles of Romberg integration can be found 
in [3]. 

8. Conclusion 

The integration method introduced in this paper provides a manner by which to approximate the definite integral 
to a high degree of accuracy without it becoming necessary to include additional information. The Romberg 
method of integration, though powerful when given a number of subintervals that is a power of two, is severely 
weakened otherwise. Once the accuracy limit is reached, the standard approach in the Romberg method is to 
double the given number of points, thereby providing an additional composite and thus an additional level of 
accuracy. 

The Romberg method of using composites that increasingly double the number of subintervals used in their 
construction, for each consecutive extrapolation, suggests that the alternative method introduced in this paper 
may be at least as accurate as Romberg's method in situations optimal for the Romberg method, but many times 
more accurate otherwise, since our alternative method provides the added flexibility of using data from other 
available composites. The Romberg method is particularly crippled when the number of subintervals available 
is odd. However, our alternative method would be successful, as composites may still be constructed in these 
circumstances. It appears that our alternative method provides a higher level of accuracy than the Romberg 
method when used on any number of subintervals that is not a power of two, and it may at least match the 
Romberg method's accuracy otherwise, and with fewer necessary approximations. 

As mentioned previously, though only four composites were used in the most thorough example in this paper, 
there is no reason to believe that this method is restricted to an error no smaller than Oiri^h^^). Providing that 
an additional unused composite can be computed, this error term should be able to be eliminated as well. The 
accuracy of this method appears to be dependent on the number of composites available. 
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